It is a testrun for v3b model checking recovery.

v3b1 model was used with standard setup (slightly wider prior on sigma + stan-specific phi calculation)

Runtime ranged 2.1-2.7k sec

3 divergences

pairs( fit, pars = c( paste0( 'r_lower[', 1:10, ']'), 'lp__'))

pairs( fit, pars = c( paste0( 'r_upper[', 21:30, ']'), 'lp__'))

pairs( fit, pars = c( 'sigma_lower', 'sigma_upper', 'lp__'))

A few divergences and funnel still clearly visible, but on the overall, it appears acceptable

pairs( fit, pars = c( 'sigma_gains', 'lp__'))

pairs( fit, pars = c( 'sigma_lpa', 'lp__'))

pairs( fit, pars = c( 'sigma_losses', 'lp__'))

pars <- c( 'sigma_lower', 'sigma_upper')
mcmc_recover_intervals( posterior[,,pars], true_vals[true_ix,pars])

Kind of overestimated standard deviations here

pairs( fit, pars = c( 'theta_upper', 'lp__'))

pairs( fit, pars = c( 'theta_lower', 'lp__'))

mcmc_parcoord( posterior[,,c( 'sigma_lower', 'sigma_upper')], np = np)

mcmc_parcoord( posterior[,,grep( 'sigma_gains', dimnames( posterior)[[3]])], np = np)

mcmc_parcoord( posterior[,,grep( 'sigma_lpa', dimnames( posterior)[[3]])], np = np)

mcmc_parcoord( posterior[,,grep( 'sigma_gains', dimnames( posterior)[[3]])], np = np)

mcmc_parcoord( posterior[,,grep( 'L_lpa', dimnames( posterior)[[3]])], np = np)

mcmc_parcoord( posterior[,,grep( 'L_gains', dimnames( posterior)[[3]])], np = np)

mcmc_parcoord( posterior[,,grep( 'L_losses', dimnames( posterior)[[3]])], np = np)

lapply( c( 'lpa', 'gains', 'losses'), function( dd) {
  apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
       c( 1, 2), function( x) {
         x %>% 
           matrix( ncol = sqrt( length( .))) %>%
           `%*%`( t( .)) %>% `[`( lower.tri( .))
       }) %>% aperm( perm = c( 2, 3, 1)) %>%
  `dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
  mcmc_parcoord( np = np)
})
## [[1]]

## 
## [[2]]

## 
## [[3]]

sapply( c( 'gains', 'lpa', 'losses'),
        function( dd) {
          apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
     c( 1, 2), function( x) {
       x %>% 
         matrix( ncol = sqrt( length( .))) %>%
         `%*%`( t( .)) %>% `[`( lower.tri( .)) %>% median()
     })
          }, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)

sapply( c( 'gains', 'lpa', 'losses'),
        function( dd) {
          apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
     c( 1, 2), function( x) {
       x %>% 
         matrix( ncol = sqrt( length( .))) %>%
         `%*%`( t( .)) %>% `[`( lower.tri( .)) %>% abs() %>% median()
     })
          }, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)

sapply( c( 'gains', 'lpa', 'losses'),
        function( dd) {
          apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
     c( 1, 2), function( x) {
       x %>% 
         matrix( ncol = sqrt( length( .))) %>%
         `%*%`( t( .)) %>% `[`( lower.tri( .)) %>% abs() %>% log() %>% sum()
     })
          }, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)

ff <- function( x) { sum( log( abs( x))) * prod( sign( x))}
sapply( c( 'gains', 'lpa', 'losses'),
        function( dd) {
          apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
     c( 1, 2), function( x) {
       x %>% 
         matrix( ncol = sqrt( length( .))) %>%
         `%*%`( t( .)) %>% `[`( lower.tri( .)) %>% ff()
     })
          }, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)

rm( ff)
sapply( c( 'u', 'l'), function( dd) {
  posterior[,,grep( paste0( '^r_', dd), dimnames( posterior)[[3]])] %>%
    apply( c( 1, 2), function( x) median( abs( x)))
}, simplify = 'array') %>%
  `dimnames<-`( c( dimnames( .)[1:2], list( parameters = c( 'upper', 'lower')))) %>%
  mcmc_parcoord( np = np)

sapply( c( 'u', 'l'), function( dd) {
  posterior[,,grep( paste0( '^r_', dd), dimnames( posterior)[[3]])] %>%
    apply( c( 1, 2), function( x) mean( abs( x)))
}, simplify = 'array') %>%
  `dimnames<-`( c( dimnames( .)[1:2], list( parameters = c( 'upper', 'lower')))) %>%
  mcmc_parcoord( np = np)

sapply( c( 'u', 'l'), function( dd) {
  posterior[,,grep( paste0( '^r_', dd), dimnames( posterior)[[3]])] %>%
    apply( c( 1, 2), function( x) min( abs( x)))
}, simplify = 'array') %>%
  `dimnames<-`( c( dimnames( .)[1:2], list( parameters = c( 'upper', 'lower')))) %>%
  mcmc_parcoord( np = np)

pars <- dimnames( posterior)[[3]]
pars <- pars[!grepl( 'pe_lpa', pars)]
mapply(function( x, a, b) (x > a) & (x < b),
       x = true_vals[true_ix,pars],
       a = summary( fit)$summary[pars,'2.5%'],
       b = summary( fit)$summary[pars,'97.5%']) %>% mean()
## [1] 0.9193254

Somewhat lower, but still acceptable recovery of parameter values (91%)